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SUMMARY 


Finite element formulations for large strain, large displacement problems 
are formulated using a kinematic description based on the corotational compo- 
nents of the velocity strain. The corotational components are defined in terms 
of a system that rotates with each element and approximates the rotation of the 
material. To account for rotations of the material relative to this element 
system, extra terms are introduced in the velocity strain equations. Although 
this formulation is incremental, in ejqplicitly integrated transient problems it 
compares very well with formulations that are not. Its simplicity, and its com- 
patibility with. constitutive equations found in "hydro" codes make it very 
attractive for this class of problems. 


INTRODUCTION 


Nonlinear structures are conventionally treated by kinematic descriptions 
that are essentially Lagrangian in nature, in that the measure of deformation 
is directly related to the total displacements. Several types of Lagrangian 
formulations are frequently used: formulations based on the Green strain or 
Almansi strain [1,2] and formulations based on corotational stretch [3,4]. 

Although velocity strain formulations have been used extensively for non- 
linear solids, as exemplified in the work of Key [5] , little study has been made 
of the application of these formulations to structures. Hughes and Liu [6] 
have presented a formulation based on the global components of the velocity 
strain. 

In this paper, a corotational velocity strain formulation will be presented 
in which the components of the velocity strain are expressed in a framework 
that rotates with the material; formulations of this type have been studied by 
Green and Naghdi [7] . The formulation is then specialized to finite elements 
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by assuming that the rotation within an element is either constant or that the 
variation in the rotation field is small or moderate. 

The potential benefits of these methods are significant. The basic equa- 
tions are simpler than Green strain or Almansi equations, which endows the re- 
sulting computer programs with both simplicity and speed. The stress conjugate 
to the corotational velocity strain is the Cauchy stress tensor expressed in 
the corotational system. Any constitutive equations based on Cauchy stress and 
velocity strain can therefore be used. Furthermore, the corotational stress 
and stress-strain matrix are both materially objective, so no Jaumann type cor- 
rections need be made for the stress state and the formulation is directly 
applicable to anisotropic materials, which is not true of the formulations 
given in [5] or [8] . 

In the next section, the fundamental equations for the corotational velo- 
city strain formulation are presented. Next, the general equations for a fi- 
nite element application of this formulation is given. In order to illustrate 
the simplicity of the method, we then give the formulation for a beam element 
assuming a constant rotation in the element. More complex relations which ac- 
count for the variation of rotation in an element are then given . The last 
section gives some examples of the application of this method to nonlinear 
transient problems. 


BASIC EQUATIONS 


We will use a kinematical and stress description by Green and Naghdi [7], 
Let us denote the material coordinates of the structure by X. , the spatial co- 
ordinates by x-j_, the displacements by u^ and the velocities fey v^. Then 


u. 
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x. 
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( 1 ) 


and the deformation gradient F^_. is given by 


F. . = 
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From the polar decomposition theorem (see [9]) it follows that the deformation 
gradient can be ejqpressed as a pure deformation, which is expressed by a sym- 
metric matrix £/^, and a rigid body rotation in the form 
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The rotation matrix FLq is orthogonal, so that 


R. . R. = 6. 

il im lm 


( 4 ) 


where 6. is the Kronecker delta, 
lm 

We will denote the coordinate system which is rotated by the rigid body 
motion of the material by x^ and call it a corotational coordinate system. 
This system is related to x^ by 


x. = R . x . (5) 

i 13 3 


and its orientation varies from point to point in the material. 

The velocity strain (rate of deformation) tensor is given by 
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( 6 ) 


and the corotational velocity strain, which is simply the same tensor with its 
components expressed in the corotational coordinates, is given by 


A 


D, 


kl 


R.. 

ik 



D. . 
13 


(7) 


or 
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( 8 ) 


The state of stress will be represented by the corotational stress T. . 
which are the corotational components of the Cauchy (physical) stress T. . 
the two are related by 
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R ik 



T. . 
13 


(9) 


265 


The corotational components of the stress are frame Invariant, so that the 
velocity strain is related to the rate of corotational stress by 


T. . 


C ijkl 



( 10 ) 


where the matrix CL . for a material depends on the state of stress and 
state variables sucfl as the yield stress but is independent of material rota- 
tion, regardless of whether the material is isotropic or anisotropic. This is 
a key advantage of corotational formulations. If the velocity strain and 
Cauchy stress are expressed in a fixed coordinate system, a Jaumann rate is 
required to provide frame invariance, but more importantly the matrix C,i j]ci' must 
also be modified to account for the rigid body rotation. Furthermore, such 
formulations are quite awkward in structural theories where it is often conv- 
venient to distinguish velocities tangent and normal to the current configu- 
ration . 

For a material in the domain fl, the rate of internal work is given by 


W = f D. . T. . dfi = f D. . T. . dfi (11) 

Jfi 13 13 ki 13 13 


finite element equations 


0 

We consider an element which currently occupies a volume . Its nodal 
displacements are nodal velocities vij and nodal forces f£j^. We represent 

the velocities within the element by shape functions 


v. 

i 


= N (x ) v 
I il 


( 12 ) 


where Nj are the shape functions which are expressed in terms of the corota- 
tional coordinates. Throughout this paper, upper case subscripts will refer to 
nodal values, as exemplified in Eq. (12) , and the indicial summation convention 
will also apply to these subscripts. 


The principle of virtual work gives 
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For elements other than the simplest, i.e. those with linear shape functions, 
the orientation of the corotational coordinates will vary within the element as 
shown in Fig. 1. Several alternatives are then available for handling the right 
hand side of Eq. (13) : 

i. a single corotational coordinate system x^ can be chosen for the ele- 
ment as shown in Fig. 1 and the relative rotations ignored; 

ii. a single corotational coordinate system x. can be chosen for the ele- 
ment and the rotations relative to x^ can be accounted for by modify- 
ing the velocity strain equations; 

iii. the relative rotations can be accounted for by using the transforma- 
tions (7) and (9) at each point of the element. 

In this paper we will explore the first and second alternative; the third 
has been explored by Hughes and Liu [6] . 

For the first alternative, the use of Eq. (12) gives 


_int 
V il f il 




T. . da 
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(14) 


so the use of the transformation 


(5) and the arbitrariness of of v^ gives 



(15) 


It should be observed that the stress is expressed in terms of a single coro- 
tational coordinate system throughout the element. Therefore, if we consider 
a beam with a constant axial stress T x , it follows that the only nonzero 
nodal forces lie along the x axis regardless of the curvature of the beam. This 
anomaly can yield spurious results whenever the flexural stiffness is small, 
since it introduces parasitic bending in states of pure membrane stress, cf 
HO] . 


The second alternative is to introduce velocity strain relations which 
account for the variation in rotations of the element but to express their compo- 
nents in the element system. If we represent these relations by 
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( 16 ) 
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then Eq. (15) becomes 


f if = \i 
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(17) 


In a subsequent section, forms of Eq. (16) for beams and axisymmetric shells 
will be presented. Higher order formulations as exemplified by Eq. (17) do 
provide better accuracy, particularly for relatively coarse meshes, but they do 
not eliminate parasitic bending. 


A SIMPLE BEAM FORMULATION 


In order to illustrate the application of a corotational velocity strain, 
we will first consider a beam element with the simplest corotational formula- 
tion where the nodal forces are evaluated by Eq. (15) . The notation used is 
shown in Fig. 2 . we will embed the element corotational coordinate within the 
element so that the x - axis always connects nodes 1 and 2. Euler -Bernoulli 
beam theory will be used, so that the velocities through the depth are com- 
pletely defined by velocities of the middle surface and V- 
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so that 


D 

x 



V y V — 
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(18b) 


where commas denote differentiation with respect to the subsequent variables. 

The velocity field V- will be approximated by linear shape functions, and 
the transverse velocity field of the midline by cubic shape functions so that 


V 

x 



(19a) 
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(19b) 


V y = £(5 3 - 2? 2 + 5) ^ + UK 3 - £ 2 ) 4> 2 

5 = x/& < 19c > 


where the rigid body part of the velocity field has been omitted since it 
causes no strain. The nodal velocities associated with deformation are thus 


{v} T = 


[V x2 ' V <f> 2 ] 


( 20 ) 


and the conjugate nodal forces 


(£ int >= tf^, V aj! 

where m^ are the nodal moments. Combining Eqs. (18) and (19), we obtain 
V A 

D x = - \ [ (6^-4) 4> 1 + (6£-2)4> 2 ] 


( 21 ) 


( 22 ) 


Equation !(15) then gives 
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where A is the cross-sectional area of the beam. 

The remaining nodal forces can be obtained from equilibrium 
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HIGHER ORDER VELOCITY 


STRAIN EXPRESSIONS 


The velocity strain expression, Eq. (22) , is exact for a beam only when the 
element's midline is coincident with the x-axis, which corresponds to the chord 
between the two nodes. If the beam element is initially curved or_if it de- 
forms enough so that the midline displaces substantially from the x-axis, these 
equations will be in error because of the following effects: 

1. The x and y axes are no longer coincident with the midline and its 

normal, respectively, so Eqs.(18) are in error because V x and V are 
not along the midline and its normal. Y 

2. The volume integration neglects the deformation of the beam relative 
to the x axis. 

In order to account for the first effect without transforming between co- 
ordinate systems within the element, higher order velocity strain relations are 
developed, using the same basic ideas employed in [4] for corotational stretch 
theories. For this purpose, the displacement of the midline relative to the 
chord, or x axis, is denoted by Y. It will be assumed that 


f = 0(e) Y,- = O(e) Y,-- = 0(e 2 ) Y,- = Y,~ (25) 


and accuracy of order e is assumed; terms of higher order and certain other 
terms are neglected. Because of space limitations, we omit the derivation, 
and give only the final result 
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For an axisymmetric shell element, we let x correspond to the radial coor- 
dinate, r, and 0 be the circumferential coordinate. The corresponding relations 
are 
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where Dg is the hoop velocity strain. The terms d^ and are identical to 
those for the beam, Eqs.(26) , while 


„ V x 1 , 

d„ = — = — (V cos a + V sm a) 

2 r r x y 


(28a) 


k = — Y, V - cos a + — [r(cos a - Y, sin a) - Y sinacosa]v - (28b) 

^ 3T X X f X IT X y j x 


where r is the current radial coordinate of a point, a the current angle between 
x and x, as shown in Fig. 2. For all applications we have considered so far, 
the second terms in the expressions for k and k 2 have been insignificant. 

EXAMPLES 


Results are given for a clamped ring shown in Fig. 3, for which experimen- 
tal results are reported in [11] ; numerical results have been reported in [3] . 
Explicit time integration with the central difference method and a lumped mass 
matrix was used. The material model is elastic-plastic with a Mises yield con- 
dition and isotropic hardening. Because the width of the ring is large compared 
to its thickness, plane strain was assumed in the z-direction. Furthermore, the 
compressibility of the elastic strains was considered neglible, so the height h 
of the cross-section was modified by 


h = h l /SL 
o o 


(29) 


for both the velocity strain computations, Eq. (18b) , and the nodal force comput- 
tations, Eq.(23) . The nodal forces were obtained by numerical quadrature using 
five points through the thickness , and two points along the length with a trape- 
zoidal rule. 

The displacements for the midpoint of the clamped ring are compared for 
plane stress and plane strain, with and without the area correction of Eq. (29) , 
in Table 1. As can be seen, the effects of the plane strain assumption are 
very significant, changing the result by 20%. The effect of the thickness cor- 
rection is less pronounced but nevertheless not neglible. 
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TABLE 1 


EFFECTS OF ASSUMPTIONS ON MIDPOINT 
DEFLECTION OF CLAMPED RING 


Assumptions 


Midpoint deflection (in) 


Plane strain, variable thickness 2.99 
Plane strain, constant thickness 3.06 
Plane stress, variable thickness 3.43 
Plane stress, constant thickness 3.53 


The plane strain solution with thickness correction compares best with 
the experiment, so we will restrict all subsequent comparisons to this case. 
The time history of the midpoint is compared to the experiment in Fig. 4. The 
reported experimental results exhibit considerable snapback, which are absent 
in the computations. Figures 5 and 6 compare the deformed shape and strain 
time histories with the experiment. Overall, the agreement is quite good and 
comparable to that of the numerical results reported in [3] . 
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DEFORMED 
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Deformed and undeformed beam element showing 
element corotational coordinates x,y and local 
corotational coordinates x,y. 
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Fig . 2 . Nomenclature 

for nodal forces 
and velocities . 
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Fig. 5. Deformed configuration of clamped ring 
compared to experimental results [11] . 
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Fig. 6. Comparison of computed strains against experi- 
mental results [11] ; angle in parenthesis is for 
the experimental record. 





